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ABSTRACT 

Aims. Radiometric ages for chondritic meteorites and their components provide information on the accretion timescale of chondrite 
parent bodies, and on cooling paths within certain areas of these bodies. However, to utilize this age information for constraining the 
internal structure, and the accretion and cooling history of the chondrite parent bodies, the empirical cooling paths obtained by dating 
chondrites must be combined with theoretical models of the thermal evolution of planetesimals. Important parameters in such thermal 
models include the initial abundances of heat-producing short-lived radionuclides (^*A1 and '"Fe), which are determined by the accre- 
tion timescale, and the terminal size, chemical composition and physical properties of the chondritic planetesimals. The major aim of 
this study is to assess the effects of sintering of initially porous material on the thermal evolution of planetesimals, and to constrain 
the values of basic parameters that determined the structure and evolution of the H chondrite parent body. 

Methods. A new code is presented for modelling the thermal evolution of ordinary chondrite parent bodies that initially are highly 
porous and undergo sintering by hot pressing as they are heated by decay of radioactive nuclei. The pressure and temperature stratifi- 
cation in the interior of the bodies is calculated by solving the equations of hydrostatic equilibrium and energy transport. The decrease 
of porosity of the granular material by hot pressing due to self-gravity is followed by solving a set of equations for the sintering of 
powder materials. For the heat conductivity of granular material we combine recently measured data for highly porous powder mate- 
rials, relevant for the surface layers of planetesimals, with data for heat conductivity of chondrite material, relevant for the strongly 
sintered material in deeper layers. 

Results. Our new model demonstrates that in initially porous planetesimals heating to central temperatures sufficient for melting 
can occur for bodies a few km in size, that is, a factor of x 10 smaller than for compact bodies. Furthermore, for high initial '"Fe 
abundances small bodies may differentiate even when they had formed as late as 3-4 Ma after CAI formation. To demonstrate the 
capability of our new model, the thermal evolution of the H chondrite parent body was reconstructed. The model starts with a porous 
body that is later compacted first by 'cold pressing' at low temperatures and then by 'hot pressing' for temperatures above x 700 K, 
i.e., the threshold temperature for sintering of silicates. The thermal model was fitted to the well constrained cooling histories of the 
two H chondrites Kernouve (H6) and Richardton (H5). The best fit is obtained for a parent body with a radius of 100 km that accreted 
at ? = 2.3 Ma after CAI formation, and had an initial '"Fe/^^Fe = 4.1 x 10"'. Burial depths of 8.3 km and 36 km for Richardton and 
Kernouve can reproduce their empirically determined cooling history. These burial depths are shallower than those derived in previ- 
ous models. This reflects the strong insulating effect of the residual powder surface layer, which is characterized by a low thermal 
conductivity. 

Key words. Solar system: formation, planetary systems: formation, planetary systems: protoplanetary disks 

1. Introduction were at work during planetary formation. Recovering this in- 
formation from these witnesses of the early history of our solar 

According to our present understanding of the formation of ter- ^y^^^^ requires to model the structure, composition, and thermal 

restrial planets, planetesimals from the size class of 100 km history of meteorite parent bodies in order to put the analytic re- 

form an important transition state during the formation of suits of laboratory investigations of individual meteorites into a 

5^ protoplanets and plan ets (e.g. IWeidenschilling & Cuzzill2006t more general context 
iNaeasawa et al.l2007h . These bodies are sufficiently big that heat 

liberated by decay of short lived radioactive nuclei, e.g. ^^Al, There are fundamental differences in the composition of or- 

does not easily flow to the surface and radiates away. Instead binary and carbonaceous chondritic meteorites that are related to 

they heat up and the bigger ones start to melt in their core region, the absence or presence of ice, respectively, in the region of the 

The pristine dust material from the accretion disc undergoes in protoplanetary disc where their parent bodies formed. We aim 

this way a number of metamorphic processes until it evolves into i" this paper to study the parent bodies of ordinai-y chondrites, 

planetary material. A few such bodies representing this interme- They were essentially ice-free, and are in this respect more sim- 

diate stage of the planetary formation process have survived in ilar to planetesimals in the inner solar system, where the terres- 

the asteroid belt and material from such bodies is available on trial planets formed. In particular we concentrate on the H and L 

earth as meteorites. Meteorites preserve in their structure and chondrites that form two rather homogeneous classes and seem 

composition rather detailed information on the processes that each to descend from a single parent body in flie asteroid belt. 

The model calculations for the thermal evolution of aster- 
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iGhosh et alJ (l2006l) . Many of the calculations are based on the 
analytic solution of the heat conduction equation with constant 
coefficients for a spherical body heated by a homogeneously dis- 
tributed and exponenti ally decaying heat source (^^Al) given by 
iMivamoto et alj (Il981h . Such models have the advantage of be- 
ing simple and easily applicable for discussions of special prob- 
lems, but they cannot be extended to more realistic cases where 
material properties like heat conduction, specific heats and oth- 
ers are neither spatially nor temporally constant, or to include 
additional physics. It was possible, however, to show that H and 
L chondrites of petrologic classes 3 to 6 originate from bod- 
ies of about 100 km size, heated by radioactive decay of short 
lived nuclei, and that the different petrologic classes correspond 
to layers at different depths within the same body which expe- 
rienced different thermal histories during the initial heating and 
subse quent cooUng of the body lO opel et al. 1994; Trieloff et all 
120031; iKleine et al.ll2008tlHarrison & Grimmll20100 . 

If more physics is to be included in a model a numeri- 
cal solution of the basic set of equations is necessary. One 
of the first models of this kind was the model calculation by 
[Yomogida& Matsui (1984) that used data for the temperature 
and porosity dependence of material properties which they de- 
termi ned by laboratory measurem ents on H and L chondrite ma- 
terial dYomogida & Matsuj|ll983h . A central point of this model 
calculation was that it addressed for the first time quantitatively 
the process of sintering of the initially strongly porous mate- 
rial under the action of pressure resulting from self-gravity once 
the body is heated to high temperatures ( > 700 K). As a result, 
the body develops a strong zoning of the material properties: A 
highly consolidated core with high thermal conductivity if tem- 
perature and pressure become sufficiently high during the course 
of the thermal evolution of the body, and a porous outer layer 
with low heat conductivity where temperatures and pressures 
never climbed to the level where compaction by sintering be- 
comes possible. 

This model shows distinct ba sic differences to a model with 
constant coefficients like that of IMivamoto et all (Il981h since 
it has an extended inner region with an only slow outwards 
drop of temperature where heat conductivity is high, and a 
rather shallow outer layer where temperature rapidly drops to 
the outer boundary temperature. The results for the depths of 
formation of different petrologic types and the predicted spa- 
tial exte nt of su ch zones ar e very d ifferent in the models of 
IMivamo to et alJ ([i981) and lYomoeTd a & Matsui (1984). This 
demonstrates the importance of including the consolidation of 
granular material by 'hot pressing' and the resulting changes of 
thermal conductivity in the models. Unfortunately the results of 
lYomogid a & Matsuil d 1 9841) are somewhat impaired by the fact 
that the important heating source ^^Al was not included in the 
calculation. In the sequel the sintering processes seem to have 
been included in the modelling of thermal evolution of a steroids 
only in the calculations of Akridge et al. ( 199% and of ISenshul 
(l2004l) . Other studies of the blanketing effect of porous outer lay- 
ers on thermal evolution ( , Akridge et al. 1998; Hey e y & Sanders; 
120061; ISahiipal et alj|2007l; iGupta & Sahiipafcoiol) assume the 
consolidation of such layers at some characteristic temperature, 
but do not include this process as part of the model calculation. 

In this paper we study the thermal history of the parent bod- 
ies of ordinary chondrites usi ng new data for the thermal con- 
ductivity of granular material dKrause et al and on the 
compaction by cold isostatic pressing before the onset of sinter- 
ing (Guttler et al. 2009). The modell ing of the sintering proc ess 
is based on the same kind of t heory dRao & Chakladei 1 1 9721) as 
in lYomogida & Matsuil dl984l) . This theory does not correspond 



Table 1. Basic mineral species considered for calculation of 
properties of chondrite material, their atomic weight A, mass- 
density g, and abbreviation of mineral name. 



species 


composition 


A 


£> 


Abbr 


Forsterite 


Mg2Si04 


140.69 


3.22 


Fo 


Fayalite 


Fe2Si04 


203.78 


4.66 


Fa 


Enstatite 


MgSiOj 


100.39 


3.20 


En 


Ferrosilite 


FeSiOi 


132.32 


3.52 


Fs 


Wollastonite 


CaSiOa 


116.16 


2.91 


Wo 


Anorthite 


CaAlzSizOg 


277.41 


2.75 


An 


Albite 


NaAlSi^Og 


263.02 


2.63 


Ab 


Orthoclase 


KAlSiiOg 


278.33 


2.55 


Or 


Iron 


Fe 


55.45 


7.81 


Irn 


Nickel 


Ni 


58.69 


8.91 


Nkl 


Troilite 


FeS 


87.91 


4.91 


Tr 



to the more recent approac hes used for modelling of hot pressing 
in technical processes (^e.gjArzt et airi 983; Fischmei ster & Arztl 
119831; iLarsson et all 119961; IStorakers et al.i il999). but it seems 
more appropriate for the lower pressure regime encountered in 
asteroids. 

The main purpose of this paper is to develop improved mod- 
els for the thermal history of parent bodies of ordinary chon- 
drites and to compare the model results with cooling histories 
of H and L chondrites of different petrologic types determined 
from thermochronological methods in cosmochemistry. This al- 
lows to better constrain the size of the parent bodies and their 
formation times. 



2. Thermal history of chondrite parent bodies 

The heat source for early differentiation and metamorphism 
of meteorite parent bod ies was decay heat of short-lived nu- 
clides like ^^Al or ^ °Fe dGonel et al.llI994i: iTrieloff et al.1 l2003t 
'Klei ne et al] 120081; iBouvier et alJ 1200 7^. If planetesimals are 
larger than tens of km, the maximum degree of internal heating is 
given by the initial ^''Al abundance, i.e., mainly formation time. 
Planetesimals formed at a; 2 Ma after calcium-aluminium rich- 
inclusions (CAIs — commonly regarded as oldest objects in the 
solar system), heat up without differentiation, yielding thermally 
metamorphosed chondritic parent bodies (ordinary chondrites, 
enstatite chondrites, or more strongly heated Acapulcoites and 
Lodranites). Primitive chondrites (such as petrologic type 3) can 
survive in the outer cool layers o f larger parent bodies (followiiig 
the onion shell model, see, e.g.. lGopel et aI.llT994t iTrieloff et al.l 
2003), o r in bodies that never grew larger than 10 - 20 km in 
size (H evev & Sandersll2006l: lYomogida & Matsuil 1 1 984 . or in 
bodies that formed relatively late. 

Such scenarios can be verified via their thermal history and 
cooling paths. These paths can be evaluated by a variety of ther- 
mochronological methods applying chronometers with different 
closure temperatures, which means the temperature where no 
parent or daughter nuclide is lost from their host minerals. The 
Hf-W dating method (e.g. Kleine et al. 2005, 2008) is useful 
for cooling below 1050 to 1 150 K. The closure tem perature for 
the U-Pb-Pb system in ph osphates is a!720 K (e.g. lGopel et al.l 
II994tlBouvier et al.ll2007h . Ar-Ar ages reflect the cooling below 
Si 550 K for oligoclase feldspar The annealing temperature for 
-'*''Pu fission tracks (e.g. Peflas et al. 1997; Trieloff et al. 2003|) 
in orthopyroxene is 550 K, for the phosphate merrillite ss 390 
K. The latter yield a relative cooling rate between 550 and 390 
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Table 2. Typical mineral composition of chondrites, mass-densities g of components, mass-fractions Xmin of min erals, and mass- 
fractions Xei of the elements that release heat by radioactive decays (data for Xmin from lYomogida & Matsuil[l983h . 



H-chondrite L-chondrite H-chondrite L-chondrite 

species composition g X^^^^ composition g X,ni„ element JCd 

g cm"^ g cm"^ 



Olivine 


Fo8oFa2o 


3.51 


Orthopyroxene 


EngjFsn 


3.25 


Clinopyroxene 


En49FSgW045 


3.09 


Plagioclase 


Ab820r6Ani2 


2.64 


Nickel-iron 


lm92Nkl8 


7.90 


Troilite 


Tr 


4.91 


Average = pbuik 




3.78 



0.37 


F075Fa25 


3.58 


0.49 


0.25 


En78FS22 


3.27 


0.23 


0.05 


En48FS8W044 


3.10 


0.06 


0.08 


Ab840r6Anio 


2.64 


0.08 


0.20 


Irn87Nkli3 


7.95 


0.09 


0.05 


Tr 


4.91 


0.05 



3.59 



Al 


9.10 X IQ-^ 


8.95 X 10" 


-3 


Fe 


2.93 X 10^' 


2.26 X 10 


-1 


K 


7.07 X 10-* 


7.07 X 10- 


-4 


Th 


5.16 X 10-** 


5.72 X 10- 


-8 


U 


2.86 X 10-** 


3.17 X 10- 


-8 



K for the respective rock (e.g. iPellas et aPllQg?!: iTrieloff et al.l 
I2OO3 V Another method to determine cooling rates between 800 
and 600 K are metallographic cooling rates, which use Fe- 
Ni diffusion profiles in metal grains consi sting of the two dif- 
ferent Fe/Ni phases kamacite and taenite dHerpfer et al.l 1 1 994t 
^opfe & Goldstein 200 L e.g.). These data are the basis for the 
simulation of the cooling of chondritic parent bodies in Sect. [8] 



3. Composition of planetesimal material 

3.1. Mineral composition of planetesimal material 

From the pristine dust material in the accretion disk, that is 
formed simultaneously with the formation of a new star, finally 
planetesimals are formed by a complicated agglomeration pro- 
cess. This process is not modelled in our calc ulation (for a re - 
view on this part of the problem see, e.g., Nagas awa et al.ll2007l) . 
We concentrate here on the evolution of the internal constitu- 
tion of the planetesimals, the successors of which can still be 
found in the asteroid belt. In particular we concentrate on the 
class of planetesimals that are the parent bodies of the H- and L- 
chondrites. Their composition is well known from studies of me- 
teorites. The matrix and chondrules that form the material of the 
ordinary chondrites consist of a mixture of minerals that them- 
selves in most cases are solid solutions of a number of compo- 
nents. We consider an average composition for the material con- 
tained in the parent bodies of the H- and of the L-chondrites that 
concentrates on the few main constituents that represent almost 
100% of the matter. The many more additional compounds with 
small abundances found in the material are not important for the 
bulk properties of those materials that determine the constitution 
of the planetesimals and their evolution. 

The pure mineral components used in our model for the 
composition of planetesimals and their mass-densities are listed 
in Table [T] The mixture of minerals that is assumed to form 
the planetesimal material is listed in Table |2l The composi- 
tion is given in the notion where, e.g., Fo8oFa2o denotes that 
80% of all basic building blocks of the mineral correspond to 
forsterite and 20% to fayalite. The quantity X^i„ denotes the 
mass fraction of the component in the mixture. The compo- 
sition of th e mixture components and the fractions Xmin are 
taken from lYomogida & Matsuil (Il983h; the average composi- 
tion of chondrites piven bv lJarosewichl(ll990h is essentially the 
same, only the Fe and Al mass-fractions are slightly less than in 
[Yomogida & Matsui ( 1983). 

The density of the solid solutions given in the table is simply 
calculated as average of the densities of the pure components. 



Table 3. Some prope rties of chon drules and matrix in H and L 
chondrites (data from lScottl (120071) ') . 







H 




L 




Chondrules 


Vol.% 


60- 


'80 


60- 


80 


Chondrule diameter (ave.) 


mm 


0.3 




0.5 




Matrix 


Vol.% 


10- 


15 


10- 


15 



weighted with their mole fractions in the solid solution. This as- 
sumes that non-ideal mixing effects are too small as that they 
need to be considered for our calculations. The average density 
g'buik of the mixture is calculated from 

^ y V' 

(1) 

V / ) 

This is the density of the consolidated meteoritic material, i.e., 
without pores. 

Table|2]presents mass-fractions X^i of the elements for which 
radioactive isotopes of sufficient abundance exist that their en- 
ergy release during decay contributes significantly to the heating 
of the planetesim als. The element abu ndances used for the cal- 
culation are from lLodders et all (l2009l) . 

3.2. Granular nature of the material 

The main constituents present in chondritic meteorites are chon- 
drules and a matrix of fine dust grains. The typical relative abun- 
dance of these constituents is shown in Table[3] they have typical 
sizes of micrometer in the case of relatively unprocessed matrix, 
and typical sizes of millimetres in the cas e of chondru les, with 
variations specific to each chondrite group ( IScottll2007 ^. Though 
chondrules may also contain fine grained mineral assemblages, 
they entered the meteorite parent body as solidified melt droplets 
after they formed by unspecified flash heating events in the solar 
nebula. 

The size of the fine dust grains of the matrix before com- 
paction and sintering of the meteoritic material is not known. 
The observed sizes of matrix particles may be not representa- 
tive because coarsening may have occurred during heating of 
the parent body (Ostwald ripening). Probably the very fine gran- 
ular units in interplanetary dust particles (IDPs) give a hint to 
the pristine size of dust particles in protoplanetary discs before 
incorporation into the forming b odies of the sol ar system. IDPs 
seem to represent cometary dust (lBradIevll20IC^ and this seems 
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to be the least processed dust m aterial handed dow n from the 
early Solar Nebula. The study of 'Rietmeijer| (1 19931) shows that 
the granular units forming IDPs are sub-micrometer sized, most 
of them having diameters smaller than 0.5 jim and many being 
nano-meter sized, <K 0.1 fi, with a small population ranging in 
size up to several jum. Typical particle radii then are roughly of 
the order of 0. 1 fim ... 0.2 fim. 

It is generally assumed that the planetesimal material ini- 
tially is very loosely packed with a high degree of porosity. As 
porosity we denote the volume fraction not filled by solid dust 
material. The volume not filled by the solids is called the pore 
space or the pores. The pores may form an interconnected net- 
work of voids and channels at high porosity, < 1, or isolated 
voids at low porosity, <«c 1 . The pores may be empty (vacuum) 
or may be filled with gas or something else. If there are no pores, 
the material is called compact. 

The compact solid material has density pbuik- This material 
is a complex mixture of minerals, that is described in Sect. 13.11 
The porous material has density 



Q - ebuik(i -4') + 



(2) 



where pp is the density of the material in the pores. If the pores 
are filled with gas then pp is small enough that we may assume 
Qp = 0. At sufficient low temperatures the pores may alterna- 
tively by filled by water and ice in which case Qj, may not be 
neglected, but here we do not consider such cases. Then 

e = ebuik(i - ■ (3) 

Besides porosity (p we also consider the filling factor 
D^l-(f>. (4) 
Then 

Q = Qhu\kD . (5) 

For this reason D is also called the relative density of the porous 
material. 

The porous material may approximately be described as a 
packing of spheres. With respect to the chondrules this is not 
unrealistic, because they show nearly spherical shape. For the 
matrix grains this may be taken as a rough approximation to de- 
scribe some basic properties of the packing. For a random pack- 
ing of equal sized spheres, Arzt ( 1982) found a relation between 
the coordination number Z, i.e., the average number of contact 
points of a particle with its neighbours, with the filling factor or 
relative density D: 



Z{D) ^Zq + C 



A) 



(6) 



Here, Zq = 7.3 and Do = 0.64. These values refer to the co- 
ordination number of a random dense packing of spheres (cf. 
lJeager&Nagelll992h . The constant C - 15.5 is the slope of 
the radial density function (distribution of centre distances). The 
approximation hold up to D > 0.9. 

For a packing of equal sized spheres it is found that there are 
two critical filling factors. One is the random close packin g with 
a filling factor of D - 0.64 , i.e., a porosity of (p - 0.36 (IScottI 
[T96llJeager & Nagell[T99l . and an average coordination num- 
ber Z - 7.3. Its formation requires taping or joggling at the ma- 
terial. The other one is the loosest close packing that is just stable 
under the application of external forces (in th e Umit of vanish- 
ing force), which has D - 0.56 ov (p - 0.44 (lOnoda & Linigerl 



Il990t LTeager & NagelllT992h and average coordination number 
Z X! 6.6. With respect to chondritic material this means that vol- 
ume filling factors of more than D - 0.64 for chondrules, like 
that given in Table |3] require that compaction by sintering or 
compaction by crushing of particles due to strong pressure oc- 
cuiTed . 

For volume filling factors of chondrules between the ran- 
dom close packing (D - 0.64) and the random loose packing 
{D - 0.56) some pre-compaction by applied forces occurred. 
The high volume filling factors of chondrules in the material 
of the chondrites, hence, shows that this material was already 
significantly compressed. It is not possible to reconstruct from 
this the initial properties of the material in the bodies that con- 
tributed to the growth of the parent bodies. We have to estimate 
this somehow. 

At the basis of the planetesimal formation process there 
stands the agglomeration of fine dust grains like that found in 
IDPs into dust aggregates of increasing size. Such an agglomer- 
ated material from very fine grains that was not su bject to any 
pressure has porosity as high as a; 0.8 . . .0.9 (e.g. lVang et al.l 
2000; Krause et al. 201 la). The number of contacts of a particle 
to neighbouring particles then is on the average equal to about 
two, i.e., the fine dust grains form essentially chain like struc- 
tures. S uch high-porosity material seems to be preserved in some 
comets (Blu m et al.ll2006l) . 

Collisions of dust aggregates during the growth process of 
planetesimals leads to comp action of the material. The experi- 
ments of Weidling et al.l ( 1200 9) have shown t hat the p orosity can 
be reduced to =s 0.64 (or even lower, see iKothe et al. 201(| 
by repeated impacts. This porosity is still higher than the ran- 
dom loose packing and, hence, is not guaranteed to be mechan- 
ically stable. The average coordination number is Z ^ 5. Lower 
porosities of <p < 0.4 were obtained by applying static pressures 
of more than lObar (Guttler et al. 2009). Since the planetesimals 
form by repeated growth process by collisions with other grow- 
ing planetesimals and impact velocities can be rather high we 
assume in the model calculations that the initial porosity of the 
material from which the parent bodies of the asteroids formed is 
already compacted to some extent and assume a porosity of the 
order of = 0.6... 0.5. 



4. Internal pressure within planetesimals 

4.1. Hydrostatic equilibrium 

The planetesimals are subject to the mutual gravitational inter- 
action between their different parts from which there results an 
inward directed gravitational force at each location. As long as 
the material is highly porous and was not yet subject to sinter- 
ing processes, the granular material may start to flow (by rolling 
and gliding of the powder particles) and evolve into a state with 
the densest packing of the granular material. If the planetesi- 
mal material is approximated by a random packing of spheres 
of equal size its p orosity in this state would be about = 0.37 
(lYang et al.l2000h . In this state the forces due to the mutual grav- 
itational attraction of all other particles are compensated by the 
reactive forces due to internal stresses that are build up within the 
particles as result of the forces acting between contact points to 
neighbouring particles. The material then can exist in a state of 
hydrostatic equilibrium with no relative motion between grains. 
This equilibrium state, however, is of a rather labile character be- 
cause at contacts points between grains stresses may be build up 
that are strong enough that particles may break apart into smaller 
pieces and some further compaction of the material is possible 
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r [AU] r [AU] 

Fig. 1. Variation of (a) mid-plane temperature and (b) pressure with distance from the star at different evolutionary epochs (0.2 Ma 
and from 0.5 Ma to 3 Ma in steps of 0.5 Ma) in a model for the evolution of the Solar Nebula (one-zone c-model). They define the 
outer boundary conditions for a planetesimal embedded in an accretion disk. In the left part the regions are indicated where the 
disappearance of an important absorber results in an extended region of nearly constant mid-plane temperature. 



by application of high pressures. For planetesimal material this 
kind of compaction is probably only relevant for impact prob- 
lems. During the gradual build-up of planetesimals, tempera- 
tures become already high enough by internal heating for sin- 
tering by internal creep to occur before really high pressures are 
built up by which the material may be crushed. 

We assume in our further considerations that during the very 
initial build-up of planetesimals the initially very loosely packed 
granular material, under the action of its own gravitational at- 
traction and/or during collisions with other planetesimals, has 
already evolved by granular flow to a state where all particles 
have a sufficient number of contacts to neighbours such that fur- 
ther relative motions are eff'ectively blocked (jp < 0.44). Further 
compaction then requires considerable pressures because parti- 
cles have to be broken for this. Additionally we neglect that the 
planetesimals may rotate and therefore neglect any deformation 
of the body resulting from this. Then we may assume that the 
planetesimals have spherically symmetric structure and the dis- 
tribution of internal pressure in the porous dust ball is described 
by the well known hydrostatic pressure equation 



dp ^ 
dr 

where 



GMy 

— ~^ 



I 

Jo 



Mr = 4n dr' r'^g(r') . 



(7) 



(8) 



The density g is the mass-density of the material. 

As long as a planetesimal is embedded in an accretion disk it 
is subject to the external pressure in the disk. The equation of the 
hydrostatic pressure stratification in the planetesimal therefore 
has to be solved with the outer boundary condition at the outer 
radius R 



p{R) = p. 



(9) 



with pext being the pressure in the accretion disk. The pressure in 
the midplane of the Solar Nebula at a number of instants as cal- 
culated from an evolution model of the disk ( Wehrstedt & Gaill 
[2002, 2008) is shown in Fig.flJ). This defines the outer boundary 
condition pext. After dissipation of the accretion disk one has to 
use 



instead. 

In order to estimate the magnitude of pressure let us consider 
a body with constant density g. We have in this case 



An -, 
Mr = —gP 



(11) 



For a body with radius R and external pressure pext integration 
of the pressure equation yields 



Pir) ^ p,.^ + ^Gg'R^[\ -[^"^ . 



P(r) = Pexl + 

and then 



For the central pressure at r = we have numerically 



PO = Pext + 1-258 X 10 



3gcm 



(12) 



(13) 



(14) 



The external pressure in the accretion disk is shown in Fig. \^ 
and we see that already for bodies as small as 1 km radius the 
central pressure in the body is significantly higher than the pres- 
sure in the accretion disk. In bodies with radii of 10 km and big- 
ger one encounters central pressures of the order of 1 bar and 
higher 

4.2. Isostatic pressing of the granular material 

The behaviour of granular material under pressure may be rather 
complex. This has been discussed by Guttler et al. (2009) with 
particular emphasis on impact processes during planetesimal 
growth. They performed also laboratory experiments for the be- 
haviour of fine grained material under the action of static pres- 
sure and how porosity is reduced if the material is compacted by 
pressing. They derived a relation between the applied pressure 
p and the relative density (filling factor) D that is observed after 
the material has come to resQ 



D(p) = 0.58 - 0.46 



1.72 



+ 1 



0.13bar. 



(15) 



p(R) = 



(10) 



Note that they denote the filling factor as ( 
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Fig. 2. Relation between relative density (filling factor) D and 
maximum experienced pressure for isostatic pressing (top) and 
run of relative density within planetesimals of the indicated size 
(bottom) that results from cold pressing due to self-gravity. 



This relation only holds for D > 0.15 because D ^ 0.15 was the 
initial value of the relative density of the material in the exper- 
iments before pressing. Its validity is also limited to pressures 
where the granular material is not yet compacted to the relative 
density D x 0.64 of a random densest packing of equal sized 
spheres. 

The meaning of Eq. ( fTSl l is that it describes the porosity 
(f) - 1 - D of a fine powder that experienced the maximum pres- 
sure p. With respect to the material in a planetesimal it gives 
the porosity distribution within the planetesimal as function of 
depth for those parts where Eq. (flST l predicts a lower porosity 
than the porosity of the surface material. The latter is determined 
by the continued impact processes during particle growth that re- 
sult also in a compaction of the powder material ( Weidli ng et al.l 
12009) . a process, however, that is not described by Eq. ( fTST l. The 
porosity <;*srf of the surface layer material has to be described in 
a different way. Hence we have to determine the porosity in the 
planetesimal from the following equation 



-I 



1 - D(p) 

<^srf 



if Dip) >l-. 
else 



'srf 



(16) 



This holds before sintering (see Sect.|6]l becomes active. Figure[2] 
shows in the upper part the variation of relative density D with 
maximum experienced pressure according to Eq. ( fTSl l. The lower 
part shows as example the distribution of relative density within 
planetesimals with 3 km to 50 km radius resulting from cold iso- 
static pressing, if no upper limit for the surface porosity is pre- 
scribed. The solutions are obtained by integration of Eq. d?) with 



Eqs. (|5]l and (flST l as equation of state. This shows that up to 
planetesimal sizes of about 7 km the material is not substantially 
compacted by self-gravity of the body. For planetesimals bigger 
than 10 km most part of the body is already compacted by cold 
isostatic pressing due to self-gravity to about a density close to 
that of the densest random packing of equal sized spheres. Only 
a rather shallow surface layer of highly porous material remains 
in such bodies. This is anyhow evident, but here a quantitative 
prediction can be made. 

4.3. Exosphere 

During the first few million years all bodies in a nascent plan- 
etary system are immersed in the protoplanetary accretion disk. 
Sufficiently massive bodies are able to gravitationally bind some 
of the gas. The minimum requirement for this is that the escape 
velocity of a gas particle from the surface of a body exceeds its 
average kinetic energy of thermal motion corresponding to the 
disk temperature 



2GMr kT 

> — 

R m„ 



(17) 



where R is the planetesimal radius and nig the average mass of 
the gas particles. From this one has for the lower limit of the 
planetesimal radius from which on gravitational bonding of an 
atmosphere starts to become possible 



R'Mm — 



IGMRin^ 
kJ 



(18) 



For a body with constant density g this means that the radius has 
to exceed a radius of 



R-, 



3kT 



^ SnGgnig 



650 



3gcm 



-3 



BOOK 



km 



(19) 



in order to gravitationally bind gas from the accretion disc and 
increase the gas pressure at their surface over the local pressure 
in the disk. This is only relevant for protoplanets with radii R > 
1000 km. 



5. Temperature structure of planetesimals 

5.1. Heat conduction equation 

For the kind of bodies that we will consider, the material has the 
structure of granular matter The internal structure of such kind 
of material is not isotropic and its properties are subject to strong 
local variations on the scale of particle sizes. As a result the 
temperature also will show local variations on the same length 
scales. The particles that form the granular material, however, 
are very small and in particular they are extremely small com- 
pared to typical length scales over which macroscopic properties 
of the bodies may vary. In this case we may average the micro- 
scopic properties of the material and also the temperature over 
volumes that contain a big number of particles and at the same 
time have dimensions small to the characteristic scale lengths 
for changes of the values of variables like average temperature 
T or average density g. Then we work only with such average 
quantities, for which we can assume that they are isotropic after 
averaging over all possible particle orientations. 

With this approximation the equation of energy conservation 
for a spherically symmetric body, expressed as an equation for 
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the temperature T of the matter, averaged over the microscopic 
fluctuations, is 



dr 1 5 2 
QCy— + — — r q,. ^ +Qh 
ar a r 



d 1 

■P—-+QV,.Fr, 

at Q 



(20) 



where c,, is the average specific heat of the granular material per 
unit mass, qr is the radial component of the average heat flow 
vector, h is the average source term for heat production or con- 
sumption per unit mass, the pdV-Xerm is the compressional work 
done per unit mass, and the last term is the work done by exter- 
nal forces. It is assumed that the external forces, F,-, are species 
independent and that no differential motion between the compo- 
nents of the material occurs (flow of gas or water through the 
porous matrix is presently not considered). The last two terms 
have not been considered so far in model calculations for thermal 
evolution of asteroids since they are negligible if the material is 
practically incompressible. However, if sintering is considered, 
the material becomes strongly compressed during the course of 
evolution and these terms have to be included. 

The quantity is the uniform radial velocity of the mixture 
components. If shrinking of the body by sintering is the only 
kind of motion of the otherwise stationary structure of the body, 
the velocity v, is obtained by differentiation of Eq. ([8]) for fixed 
Mr with respect to time. There foUows 



-2 f^^'^' 

Jo 



2 ^ 
dt 



(21) 



We will consider models of planetesimals that are in hydro- 
static equilibrium without internal flows. There may be some 
very slow radial motion of the matter if the material starts to 
shrink by sintering at high temperatures. This kind of extremely 
slow motion is completely negligible in the substantial deriva- 
tives d/df - d/dt+Vrd/dr in Eq. ( |20] |. However, one consequence 
of shrinking is not negligible: If the body shrinks the gravita- 
tional potential energy decreases and the corresponding amount 
of energy is transferred to the matter as heat. This is described 
by the term gVrFr where Fr is the local gravitational acceleration 



g = — 



(22) 



The term corresponding to the work done by the forces has to be 
retained. With this approximation we have the following equa- 
tion for the temperature 



dT 1 d 



d 1 



QCv -7— + -J — r qr ^ +gh - P - QVr 



GMr 



dr 



f2 Qf 



dt g 



(23) 



The variation of q with time is discussed in Sect. |6l 

The heat flow vector has contributions from a number of pro- 
cesses. For the solid component of the material there is a con- 
tribution from heat conduction by phonons or, in the case of 
electric conductors (e.g. iron), from conduction electrons. For 
a porous material there is also a contribution from the heat con- 
duction by the gas in the pores. If the material is translucent then 
one has also to consider a contribution to heat conduction by ra- 
diative transfer. All these processes have the property, that their 
contribution to the total heat flow is proportional to the gradient 
of the temperature. Generally the coefficient of proportionality is 
a second rank tensor, except if the properties of the material are 
isotropic, in which case it degenerates to a simple scalar factor 
For granular material the local transport properties for heat are 
by no means isotropic. We assume, however, that after averag- 
ing the average heat flow vector is proportional to the gradient of 



1200 



1000 




1600 



Fig. 3. Specific heat of planetesimal material for mineral mix- 
tures of H chondr ites and L c hondrites. Dotted line is the analytic 
approximation of lYomogida & Matsui ( 1984). 

the average temperature. The radial component of the heat flow 
vector then takes the specific form 



dT 

qr^-K — 

or 



(24) 



with some average heat conduction coefficient K that is different 
for each of the different transport processes. 

The essential material properties that enter into equation ( l23l l 
are the specific heat per unit mass, c,,, the heat conductivity K, 
and the heat production. In the following subsections we de- 
scribe how we determine these quantities for the material of the 
parent bodies of ordinary chondrites. 

The body is also subject to heat and matter exchange with its 
environment. This is treated by defining appropriate boundary 
conditions for Eq. (l23T l that are discussed after our discussion of 
the material properties. 

5.2. Heat capacity of material 

The heat capacity c,, of a mineral mixture is simply the weighted 
sum of the heat capacities of its components. It is calculated in 
our model calculations from 



(25) 



where Xmin,, is the mass-fraction of the /-th component in the 
mixture of solid components and c,,., the heat capacity per unit 
mass of component /. The quantities Xmin.; are given in Table |2] 
Since we intend to consider bodies of a size of no more than a 
few 100 km, pressures remain below the kbar-range (see Eq.[T4li 
and compression of solid material under pressure is negligible 
because of the low compressibility of minerals. Under these con- 
ditions there is no significant difference between Cp and c,, and 
we may use Cp instead of c,,. For our model calculations data for 
Cpj(T) are taken from the compilation of Barin ( 1995) that gives 
the heat capacities Cp per mole. These quantities are converted 
to heat capacity per unit mass by dividing by the mole mass M,- 
of species /: 



(26) 
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Table 4. Data for calculating heating rates by decay of radioac- 
tive nuclei. 



10^ 



Species 


/ 


E 

[iVlc V J 




T 

[al 


h 

W Kg 




2«A1 


5.1 X 10-^ 


3.188 


(1) 


1.0 X 10« 


1.67 X 10- 


7 


sope (high) 


1.6 X 10-"^ 


2.894 


(1) 


3.8 X 10'' 


2.74 X 10- 


8 


(low) 


4.2 X 10-' 












40k 


1.5 X 10-3 


0.693 


(2) 


1.8 X 10' 


2.26 X 10- 


11 


232-ph 


1.00 


40.4 


(2) 


2.0 X 10'° 


1.30 X 10- 


12 


235u 


0.24 


44.4 


(2) 


1.0 X 10' 


3.66 X 10- 


12 


238u 


0.76 


47.5 


(2) 


6.5 X 10' 


1.92 X 10 


12 



Sources: (D IFinoccM &"GaiilfT99l . (2^ IVan Schmu3 fT99l 

Values of c,, for the required temperatures are determined by in- 
terpolation in the tables for Cpj. The heat capacity for solid solu- 
tions is calculated as weighted mean of the heat capacities of the 
pure components taking their mole fractions as weights. 

The variation of the specific heat of the mixture is shown in 
Fig. [3] Since some of the minerals suffer structural transitions at 
certain temperatures and since Cp of iron has a cusp at the Curie- 
temperature of 1042K, the te mperature variation of Cp shows 
some kinks and jumps (cf. also Ghosh & McSween 1999). They 
might be sources of numerical problems. For comparison the fig- 
ure also shows th e analytic approximation fo r Cp for the bulk 
material given bv 'Yomogida & Matsuil (11984'). This is an alter- 
native if the jumps in the temperature variation of Cp result in 
numerical problems, but in our calculations it was not necessary 
to take recourse to this approximation. 

5.3. Heating by radioactive nuciei 

Next we consider the source term h in Eq. ( |23] ). There are es- 
sentially two different kinds of sources and sinks of heat within 
the planetesimal bodies. One source is the energy liberated dur- 
ing decay of radioactive isotopes of a number of elements. The 
other one is consumption of energy during melting of planetes- 
imal material or liberation of energy during solidification of the 
melt. Melting is not considered because we aim to study parent 
bodies of undifferentiated meteorites. 

The main sources of heat input by radioactives during the 
early heating up of planetesimals and the subsequent cooling 
phase are decay of ^^Al and possibly ^"Fe (cf. the discussion in 
iGhosh et al.l l2006l) . More long-lived radioactive nuclei, essen- 
tially -^^^Th, 235,238y^ ^jj^ 40j^^ ^j.g responsible for the later long- 
term evolution of the temperature. For the nuclei decaying by 
jS-decay (^^Al, *"Fe, "^''K) the energy of the fast electrons and of 
the emitted y-photons is absorbed within the planetesimal mate- 
rial and is converted to heat, while the neutrinos leave the bodies 
and carry away their part of the energy. For the nuclei decaying 
by a-decay the whole decay energy is absorbed by the planetesi- 
mal material. We assume that no chemical differentiation occurs 
in the bodies that we consider. Hence, after having averaged over 
the inhomogeneous microstructure of the material, the heat pro- 
ducing nuclei are homogeneously distributed over the bodies. 
The heat production rate by these nuclear processes is 



^nuc — / , fi e ^ 



(27) 




100 



t[Ma\ 



The sum is over all nuclei that contribute to heating, Xei_, denotes 
the mass-fraction of the corresponding element in the material of 



Fig. 4. Contributions of different heating mechanisms to the to- 
tal heating rate. Time f is after formation of CAIs. The dashed 
line indicates the formation time of the H chondrite parent body 
(2.3 Ma, see Sect.|8]l. The release of gravitational energy by con- 
traction is shown for the final model of the H chondrite parent 
body. 

the planetesimals (see Table IS, mei,; the atomic mass of the ele- 
ment for the isotopic mixture at the time of formation of the so- 
lar system, f, the fraction of the isotope of interest at the time of 
formation of the solar system, r,- the decay time scale for e-fold 
decrease of abundance of the isotope, and t is the time elapsed 
since formation of the solar system. 

The constants for calculating /znuc are given in Table |4] The 
element abundances used for calculating Xei., for the minera l 
mixture given in Table |2] are taken from Lodd ers et a n (I2009h . 
Isotopic abundance s of K and U at time of solar system forma- 
tion are taken from (Anders & Grevessd (fT989). The abundance 
of ^''Al is that given by Nvauist et al. (2009). The abundance of 
^"Fe is disputed in the literature. Table |4] gives the probably up- 
permost value from Birck & Lugmair ( 1988) and the lower limit 
according to Ouitte et al. (2007). There are indications that ^°Fe 
was not hom ogeneously distributed in the early solar system 
(lOuitte et al.l l2010). which means that the initial ^°Fe abundance 
in the parent bodies is not known a priori and is an additional 
free parameter for the modelling. For the decay time of ''"Fe the 
recent revised val ue for the half-lfe ti/2 = 2.62 + 0.04 Ma of 
'Rugel et all ilQQ^ is used. 

Figure|4]shows the variation of the heating rate per unit mass 
with time elapsed since formation of CAI, calculated with the 
high ''"Fe abundance (see Table|4]for this). The dominant heating 
source is ^*'A1 at the time of formation of planetesimals (f < 
5 Ma), but ''"Fe dominates as heat source for an extended period 
from a: 5 Ma to a: 20 Ma because of the revised ^"Fe half-life. 

For comparison Fig. |4] also shows the contribution of the re- 
lease of gravitational energy to heating during shrinking of the 
body, resulting for the model of Sect. [8] For the rather small 
bodies that are considered in this paper this heating source is not 
important (but included in the model). 

5.4. Heat conduction by tiie porous solid material 

For the heat conductivity K of the chondritic material we use two 
different types of experimental data. For low porosities from the 
range < <p < 0.2 we use d ata measured for a number o f or- 
dinary H and L chondrites bv lYomogida & Matsuil (Il983l) . For 
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Fig. 5. Variation of heat conductivity K with porosity 4>- Results 
for fine grained silica powder (filled circles) from experiments of 
Krause et al.' ("201 la'.'W), and for particulate basalt (crosses) from 
Fountain & West ( 1970). Typical grain size is indicated for both 
cases. Open squares and open circles are experimental results 
for heat conductivity and porosity for ordinary H a nd L chon- 
drites, respectively, from dYomogida & Matsuill l983h . Solid fine iKrause et~ar 
is analytic fit, Eq. ( l30l l. to the data. 



1201 Ibi slightly different values of the constants were given). 
This fit is shown as the second dashed line in Fig. |5] 

The experiments of Krause et al. (201 la) are conducted with 
very fine grained silica powder This is not the same kind of ma- 
terial as it is found in chondrites, but for two reasons it may be 
considered as a reasonable proxy for chondrite material before 
strong compaction. First the mineral mix in chondrites is domi- 
nated by silicates and all silicates have similar heat conduction 
coeflicients. Second the heat conduction of very loosely packed 
material is via the tiny contact regions between the particles. In 
chondrites part of the material are the rather big chondrules (size 
a: 1 mm), but as long as the material is not strongly compacted 
and the chondrules are well separated by the very small grained 
matrix (particle sizes < 0.5 /im), the heat conduction is obvi- 
ously governed by the heat flow through the contact points be- 
tween the tiny matrix particles. In this respect the basic physics 
of the heat transport in the experiment of Kra use et al.l ( 1201 lal) 
should be very similar to that in chondrite material before strong 
compaction. 

The powder particles used in the experiments are about a 
factor of five to ten times bigger than the matrix particles in 
chondrites (iRietmeiieill 19931) . Some indications on the influence 
of particle sizes may be obtained by comparing the results of 



zes may 
(1201 lah ^ 



^ with result s of heat conduction measure- 

ments of )Fountain & West! (Il970h for powders of basaltic partic- 



high porosities (p > 0.4 we use the data for silica powder de - 
rived from laboratory measurements by iKrause et al. I (l20TTi) . 
All these measurements were conducted under vacuum condi- 
tions in order to exclude any contribution from heat transport 
from gas-fillings in the pores. 

Figure |5] shows conductivity K plotted versus porosity for 
chondritic material at T = 300 K. The data for H and L chon- 
drites scatter significantly and because of the small number of 
available data points no obvious systematic difference between 
the two types of material can be recognized. Therefore we fit 
both sets of data with a single analytic approximation 



(28) 



with two constants K\, and (pi. This exponential form enables 
a reasonable fit of the data. The constant K\, may be inter- 
preted as the extrapolated average thermal conductivity of the 
bulk material at vanishing porosity, for which we obtain - 
4.3 Wm~ ' K~'. For the seco nd constant we choose <pi = 0.08 
(see also IKrause et al. l l20Tl b). In Fig. (H the data are normal- 
ized with the value of K\,. The dashed line running through the 
data points for the chondrites give our approximation Ki((p). 

At high porosities Fig. (|5]l shows the conductivity K for a 
silica powder that consists of equal sized spheres with 1.5 /im 
diameter. By the n ature of the experimental method the data of 
IKrause et al.l (|201 la) do not refer to a well defined temperature 
but the heat conductivity was derived by analysing the cooling 
behaviour of their sample. Therefore the value of K is some aver- 
age value over the raise and fall of the temperature in the exper- 
iments, which is somewhat above room temperature. The data 
are fitted with an analytic approximation of the form 



ulates that are much coarser grained. Figure |5] shows results for 
their size separate with average particle size of 50 um. Th e vari- 
ation of K with porosity for the iFountain & West! d 19701) gran- 
ular material is very similar as of the silica powder used by 
IKrause et al.l (l2011ah . except that the conductivity is lower by 
a factor of somewhat less then a factor of ten. The particle sizes, 
on the other hand, are bigger on average by more than a factor 
of t hirty. Thi s suggests that the heat conductivity measured by 
Kra use et al.l |201 la) for the highly porous silica powder under- 
estimates the conductivity of the matrix material of chondrites at 
the same porosity, but probably not by a big factor. In case of a 
power law dependence of K on particle sizes one may speculate 
that a by a factor of about three higher conductivity of meteoritic 
material than for the silica powder may be an appropriate esti- 
mation, but with out mo re definite information we take the values 
of Krause et a fl (l20TTi) for our model calculations. 

The two fits, Eq. (|28]l and Eq. (|29]l, for < 0.2 and > 0.4, 
respectively, are combined into a single analytic approximation 
for K{^) by 



1/4 



(30) 



a-0/02 



(29) 



with two constants a and (pi and the same value of Ki, as before. 
This type of exponential dependence on (p allows a reasonable 
fit of the data points also in this case. The constants are found to 
be a = 1.2 and 4>2 - 0.167 (in an earlier version. iKrause et"an 



in order to smoothly interpolate between the two limit cases, in 
particular in the intermediate range of porosities 0.2 < (p < 0.4. 
This approximation is shown as the full line in Fig.|5] 

The bulk conductivity A^b in Eq. (l28T l is temperature depen- 
dent. Fig ure [6| shows data fo r K fo r H and L chondrites as 
given by 'Yomogida & Matsuil (Il983l) . divided by Ki ((p), given 
by Eq. (l28l l. The data for K{T) / Ki{(p) show for each of the mete- 
orites a clear systematic variation with temperature. The extent 
of these variations, however, is much less than the scattering be- 
tween the different meteorites which amounts to variations by a 
factor of about two. The origin of these meteorite-to-meteorite 
variations is not known, but most likely they have their origin 
in variations in the composition and structure of the meteoritic 
material. Such variations can presently not accurately be ac- 
counted for and will not be considered in our model calculations. 
Therefore we will also neglect the small variation of K with T 
and take Ki, in our calculations to be temperature independent. 
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Fig. 6. Temperature variation of thermal conductivity K. Plotted 
are data for H chondrites (open rectangles) and L chondrites 
(open circles), normalized with the ^-variation according to ap- 
proximation (l28b . The lines connect data for each of the mete- 
orites. 



5.5. Heat conduction by other processes 

The material of the planetesimals is dominated by minerals that 
are transparent in the far infrared region. Therefore energy trans- 
port by radiation is possible. The energy flux by radiation in- 
creases strong with increasing temperature and for temperatures 
of the order of about 800 K and higher the heat flux by radiation 
becomes non-negligible. The complicated structure of the mate- 
rial (chondrules densely packed in a porous dust matrix) makes 
it difficult to calculate this from first principles. We follow there- 
fore the proposal of Yomogida & Matsui ( 1984) to take the re- 
sults of the laboratory measurements of Fountain & Wes% (1 19701) 
of radiative heat conduction in basaltic powders as an approxi- 
mation for meteoritic material. The model calculations show that 
for the bodies of interest (temperature stays below melting tem- 
perature) radiative heat conduction never becomes an important 
energy transport mode. This is because if temperature becomes 
high enough for radiative heat conduction to contribute some- 
what to the net heat flux, the onset of sintering above s; 700 K 
results in a strong increase of heat conductivity by phonons that 
outnumbers radiative contributions again. 

A possible contribution to heat conduction by pore gas is 
small and neglected in this paper; details of this will be discussed 
elsewhere. 



5.6. Boundary conditions 

For solving the heat conduction equation ( l23T l one has to specify 
an initial condition T{r, to) at some initial instant fo and boundary 
conditions, in our case at r = and r - R. 

The inner boundary condition is that there is no point source 
for heat at the centre which translates into the Neumann bound- 
ary condition 



dT 
57 



0. 



(31) 



r=0 



The temperature of the surface layer of the body is deter- 
mined by the equilibrium between (i) the energy fluxes toward 
the surface from the interior and from the outside, and (ii) the 
energy fluxes away from the surface to exterior space (cf., e.g.. 



iGrimm & McSweenI [19891 iGhosh et al.ll2003h . There results a 
mixed boundary condition at the surface 



-K 



dT 
57 



(32) 



The left hand side describes the energy flow from the interior 
toward the surface by heat conduction. The first term on the rh.s. 
is the rate of energy loss by radiation from the surface to exterior 
space, the second term, Fgxt, is the rate of energy gain by outer 
sources. During the first few million years, if the planetesimals 
are still embedded in an optically thick accretion disc, Fext is 
given by osbT'^, with being the temperature at the midplane 
of the disk (see Fig. [T^). After disk dispersal the planetesimal is 
irradiated by the proto-sun and the rate of energy input, Fext, is 
given by (l-Asuif)crsBr^^J/4fl^. Here a is the (average) distance 
to the star and Asmf is the albedo of the planetesimals surface. 

Alternatively one can consider at the surface the Dirichlet 
condition 



T{R) = Tb 



(33) 



with some prescribed value T\,. This has been done in a number 
of published model calculations where some fixed value for 
was assumed. 



5. 7. Initiai condition 

Large planetesimals with radii of the order of 100 km occur 
as transition states during the growth from km sized bodies of 
the initial planetesimal swarm to protoplanets with sizes of the 
order of 1 000 km. The growth initially proceeds rather slowly 
on timescales of a few 10^ years until run-away growth com- 
mences after the bigg est planetesimals reach sizes of about 
10-20 km (e.g. iWeidenschilling & Cuzzill2006t iNagasawa et al.l 
l2007h . During run-away growth the mass rapidly increases 
within less than 10^ years to the size of a protoplanet. This means 
that the bodies are formed on timescales shorter than the decay 
time of 1 Ma for ^^Al and also that they collect most of their 
mass within a period even much shorter than this. There is not 
sufficient time available for strong heating by ^^Al decay of those 
planetesimals that contributed to the growth of a 100 km sized 
body; most of the heat released by radioactive decay is released 
after its formation. 

Therefore we will base our calculations in this paper on the 
"instantaneous formation" approximation where it is assumed 
that the body is formed within such a short period that all heat- 
ing occurred after its formation. Within the framework of the 
instantaneous formation approximation it is appropriate to pre- 
scribe the temperature Tc of the disk material at the formation 
time of the planetesimals as initial value of the temperature. 

Modifications resulting from a finite duration of the growth 
period will be considered elsewhere. 



6. Sintering 

The compaction of material in planetesimals is a two-step pro- 
cess. The initially very loosely packed dust material in the plan- 
etesimals comes under increasing pressure by the growing self- 
gravity of the bodies. The granular material can adjust by mu- 
tual gliding and rolling of the granular components to the ex- 
erted force and evolves into configurations with closer packing. 
The ongoing collisions with other bodies during the growth pro- 
cess enhances this kind of compaction of the material. This mode 
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of compaction, "cold pressing", by its very nature does not de- 
pend on temperature and operates aheady at low temperature; 
it was considered in Sect. 14.21 A second mode of compaction 
commences if radioactive decays heat the planetesimal material 
to such an extent, that creep processes are thermally activated in 
the lattice of the solid material. The granular components then 
are plastically deformed under pressure and voids are gradually 
closed. This kind of compaction by "hot pressing" or "sintering" 
is what obviously operated in ordinary chondrite material and 
the different petrologic types 4 to 6 of chondrites are obviously 
dif ferent stages of compaction by hot pressing. 

[Yomogida & Matsui (1984) were the first to perform a quan- 
titative study of this process by applying early theories of sin- 
tering developed in material science. We follow here essentially 
the same approach because more advanced modern theories of 
hot pressing are developed to model metallurgical processes that 
apply generally much higher pressures (» 1 kbar) than what 
is typically encountered in compaction of planetesimal material 
(< 10^ bar) and are mainly concerned with the final stages of the 
process. Because the rate of increase of temperature in planetesi- 
mals is very low (of the order of 1 0"^ K yr" ' ) the creep processes 
result in finite deformations of the material already at rather low 
temperatures or pressures, where under laboratory conditions no 
effect is seen. The more simple early theories of hot pressing 
seem better to fit to such situations. 

6.1. Equations for hot pressing 

For describing the sintering process we initially assume a dense 
packing of equal sized spheres with initial radius Rq. The pack- 
ing is sufficiently dense that no further compaction can be 
achieved by pressing without crushing the spheres. On average, 
the individual granular units will touch each other at Z contact 
points. At sufficiently high pressure and temperature the individ- 
ual spheres will plastically deform at the contact points by creep 
processes and contact faces will form between adjacent particles 
while the volume of the particle will remain constant. As sinter- 
ing proceeds, the voids between the spheres become smaller and 
the sphere centres get closer 

There are two stages for this process. In the first stage the 
voids form an interconnected network between the granular 
units. This closes at some stage of the sintering process and there 
remain isolated pores, that have to close by further sintering in 
a second stage. The relative density D at the transition between 
stages one and two depend on the type of packing. The following 
approximations are for the first stage. 

The first basic assumption of the de formation theory of hot- 
pressing bxJKakar & Chakladd dUl^), on which the work of 
lYomopda & Matsui (1984) is based, is a pure geometrical one. 
It is assumed that the formation of the contact faces can be con- 
ceived as if at each contact point a cap would have been cut-off 
from each of the two contacting grains. Then, for grains of equal 
radius, the contact areas are circular areas with radius a. It is as- 
sumed that all cut-away caps have the same height h and radius 
a at their base. The volume of one such cap is 



o 

and its height 



h^R- V/?2 - fl2 . 



(34) 



(35) 



The granular units then are (by assumption) spheres with Z caps 
cut-off from them. To conserve the original volume of the sphere. 



the radius R of such a truncated sphere has to be bigger than the 
pristine radius Rq. Conservation of volume requires 



4;r 3 47r 3 



(36) 



This holds as long as the contact areas do not come into con- 
tact with each other. The relation to the relative density D is 
(iKakar & Chakladed[T967h 



(37) 



Here Do is the relative density of the initial packing of spheres 
with radius Rq. For a given number of contact points Z and given 
D{), Rq, Eqs. ( l34l l to dSTl i define R and a in terms of the relative 
density D. 

In the theory of iRao & ChakladeJ (Il972h a number of reg- 
ular packings of spheres is considered for which the number of 
contact points Z is fixed (Z - 6, 8, 12). In particular they favour 
the 'hexagonal prismatic' packing with Z = 8 and gave their 
formula for this case. This is the model that has been used by 
Yomogi da & Matsuil d 19841) in their study of sintering of plan- 
etesimals. They argued that many experiments on the packing of 
small spherical particles of constant size show that the porosity 
achieved after sufficient tapping would be near 40%, with an av- 
erage of about eight contact points per particle. Since a regular 
hexagonal prismatic packing of spheres also has a coordination 
number of eight and a porosity of 39.5%, and a random close 
packing has porosity 36% and on average Z = 7.3 (see Sect. 13. 2b . 
they used that packing model in their sintering models. For a dis- 
cussion of t he more ge neral ca se of random packing s see, e.g., 
lArzi (ll982h : lArzt et all (ll983h : lFischmeister & Arztl(ll983h : the 
equations obtained in that case are more involved, but there are 
no basic differences. 

The second basic assumption in the theory of 
iRao & Chakladej (Il972h of hot pressing is that the strain 
rate is related to the applied stress by the relation for power-law 
creep 



Act?, 



(38) 



and that e is given in terms of the rate of change of relative den- 
sity as 



D 

"d 



(39) 



The stress cri is the pressure acting at the contact faces of the 
granular units. The quantities A and n have to be determined ex- 
perimentally for each material. The quantity A depends on tem- 
perature. 

The stress ct] is given by the pressure acting at the contact 
areas between the granular units. It is assumed that this is given 
in terms of the applied pressure p and the areas of contact faces, 
na^, and average cross-section of the cell occupied by one gran- 
ular unit, Cav, as 



(40) 



In the hexagonal prismatic packing model favoured by 
lYomogida & Matsuil (Il984l) . the cross-section Cav is given by 



Cav =2V3(7?2-fl2)_ 



(41) 



Via the dependence of R and a on D this is a function of D . 
Values for other packing models (e.g. Ka kar & Chakladeill 1967b 
are within 0(1) of this. In the initial stages of sintering (small 
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a) the stress cri is much higher than p, in the final stages it ap- 
proaches p. 

Equations (l38l l. (39[ . (l40t result in the following differential 
equation for the relative density 



dD ^./Cav V 

— ^ -DA[—p\ 



(42) 



With the above relations between R, a and D this is a (closed) 
set of equations for calculating the time evolution of D, that has 
to be solved together with the other equations for the structure 
and evolution of the planetesimal which define the pressure and 
temperature. 

The transition to stage 2 by closure of voids is a ssumed in 
models of hot pressing to occur at D > 0.9 (e.g. lArzt et al] 
1983h. The equation for D becomes oc 1 - D for this case (cf. 
Wilkinson & Ashbvlll975t lArzt et al.l[T983h . Because the corre- 
sponding full equations are similar in structure to -Eq. ( l42l i ex- 
cept for the factor 1 - D we include the final pore-closing stage 
simply by multiplying for D > 0.9 the r.h.s. of Eq. (l42T i by a 
factor 

F=10-(l-D) (D>0.9), (43) 
to get a continuous transition between both cases. 

6.2. Data for olivine 

The pre-factor A and the power n in Eq. ( 138)1 have to be de- 
termined by la boratorv exponents. .Yomo gid a^~Matsuil (1 19841) 
used data from lSchwenn & Goetzd (Il978l) for olivine. No other 
data for materials of interest for planetesimals seem to have been 
determined since then. Schwenn & Goetze (1978) gave the fol- 
lowing fit to their experimental data for small spheres of olivine 
{R < 53, //m): 



«3 



(44) 



with cr stress on contact faces (in bar), E-^^i the activation energy 
for creep, T the temperature, /?gas the gas constant, and R the 
radius of granular units (in units cm). 

For the ac tivation energy a value of iS act - 85 + 29 kcal mol"^ 
was given bv ISchwenn & Goetzd (11978'). In the model calcula- 
tions we use the value Eact = 85kcalmol For the pre-factor 
A a range of values from 1 .6 x 10"^ to 5.4 x 10"^ was given by 
ISchwenn & Goetz"3 ( Il978h : as a compromise we use a value of 
A = 4 X 10~^ in our model calcula t ions. 

Note that Yomogida & Matsui (1984) choo se to use for e the 
approx imation with « = 1 given in Eq. (7) in S chwenn & Goetzd 
dMi, wh ile we prefer to use use the approximation given by 
Eq. (8) of ISchwenn & Goetzd (Il978l) . because they explicitly 
state that this describes their measured cr-dependence. 



7. Results for thermal evolution of planetesimals 

7.1. Model calculation 

The calculation of a model requires to solve the differential 
equations for heat conduction, Eq. ( l23T l, hydrostatic equilib- 
rium, Eq. (|7]), for the evolution of porosity, Eq. ( l44l l. together 
with equations for the material properties, the equation of state, 
Eq. (|3]l, the equations for heat conductivity, Eq. (l30b . and heat ca- 
pacity, Eq. ( 1261 ), and together with appropriate initial and bound- 
ary conditions. 



Table 5. Model parameters for the model of iMivamoto et all 
(1198 1*) (MFT81) for a consolidated body (average L chondrite), 
and for a similar model of an initially porous body without (PLO) 
and with (PL) additional heating by decay of ^"Fe and long-lived 
nuclei. 



Quantity 




MFT81 


PLO 


PL 


Units 


radius 


R 


85 


85 


85 


km 


formation time 


^form 


2.4 


2.3 


2.3 


Ma 


heat conductivity 


Kb 


1 


4.3 


4.3 


Wm' K' 


surface temperature 


T, 


180 


150 


150 


K 


density 


Qbu\k 


3.2 


3.59 


3.59 




-«A1/2'A1 




5 


5.1 


5.1 


xlO^^ 


'^"Fe/^'^Fe 









4.1 


xlO-^ 


initial porosity 




(10%) 


60% 


60% 





The heat conduction equation and the pressure equation are 
re-written in terms of Mr, defined by Eq. as independent 
variable and are discretised for a set of fixed mass shells with 
masses AM, (/ = 1, /) and shell boundaries r, (/ - 0, /). The Mr- 
coordinate corresponds to a Lagrangean coordinate that is fixed 
to the matter. For this choice of coordinates there is no flow of 
matter across cell boundaries. This enables a simple treatment of 
growth of the body, if this is considered, and it avoids problems 
with numerical diffusion in case of inhomogeneous composition 
(e.g., radial variation of porosity). 

The heat conduction equation is solved by a fully implicit 
finite difference method with Neumann boundary condition, 
Eq. dSTJ, at centre and Dirichlet boundary condition, Eq. ( 1331 ). at 
the surface. The first order ordinary differential equation for 0, 
Eq. (l44b . is solved by the fully implicit Euler method. 

In order to account for the non-linear coupling between the 
different equations we perform a fixed-point iteration where we 
solve the equations in turn as follows: 

1 . Given are values of <pi, Tj, AM, for each mass shell / at some 
instant t^-i. We have to calculate new values at next instant 
tk - tk-i + Af. The values of 0,-, T/ at t^-i are used as starting 
values for the iterative calculation of <pi, T,- at fj. 

2. The heat production by radioactive decays over the period At 
is calculated for each shell. 

3. For given porosity 0,- one finds qi from Eq. (O and with given 
mass AM,- in each mass shell / we calculate, starting from the 
centre, the shell boundaries r, at f^^. 

4. From the change of r, over time At we find the grid velocity 
and the heat production by release of gravitational energy for 
each shell (last term in Eq.l23Tl. 

5. We calculate the pressures p, at shell boundaries from the 
discretised pressure equation, starting with the given external 
pressure at the surface. 

6. We solve for given temperatures T, and pressure pi at each 
grid-point Eq. (l42l i for the porosity over time interval Af to 
determine an updated value of 0, at f^. The corresponding 
non-linear equations are solved iteratively with an accuracy 
of better than 10"'^. 

7. We calculate from the updated porosity and pressure the heat 
conductivity. 

8. We calculate for given T, the heat capacity. 

9. The surface temperature is determined from Eq. (|32] |. This 
equation is solved for T^, using on the l.h.s. the values for T, 
and K from the last iteration step. 
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Fig. 7. Temperature evolution of a body of 85 km radius at dif- 
ferent depths from the surface and at centre: Dotted line at depth 
4.25 km, short dashed line at 17 km depth, long dashed line at 
23 km depth. Full line shows temperature evolution at the cen- 
tre, (a) The model is cal culated with the same p hysical input as 
in the analytical model of'Miva moto et al.l(ll981 b. Model param- 
eters MFT81 in Table |5] (b) Similar model, but now calculated 
for a porous body, considering thermal conductivity of porous 
material according to Eq. ( l30l l, and sintering and cold pressing 
as described in this paper. Model parameters PLO in Table |5] 
(c) Same kind of model as (b), but now additional heating by 
^°Fe and long-lived nuclei considered. Model parameters PL in 
Table S 



10. Updated values of temperatures T, at are calculated for 
all / from the difference equations resulting from the heat 
conduction equation. 

11. Check, if deviation of new values of T, and 0, at from 
current values is sufficiently small. 

12. If not, repeat calculation from step 3 on. 

13. Otherwise determine new stepwidth Af (see below), advance 
k by one, and repeat from step 1 on for the next time step. 



This kind of iteration converges usually within about ten to 
twenty steps. The accuracy requirement was that the relative 
change of T,, <pi between subsequent iteration steps is less than 
10 This simple scheme works efficiently since the coefficients 
in the heat conduction equation do not strongly depend on tem- 
perature. Then this equation can be solved by the simple fixed- 
point iteration described. Test calculations done with complete 
linearisation of the non-linear equation showed that this did not 
significantly improved the efficiency of the solution method in 
our particular case. The advantage of our method is that it poses 
less stringent requirements on the existence and continuity of 
derivatives than the Newton-Raphson method for convergence 
of the iteration method. 

The boundary condition given by Eq. ( l32l l in principle should 
be built into the difference equations for the heat conduction 
equation. Numerical experience showed that this occasionally 
resulted in stability problems. Our present method is to solve 
Eq. ( |32] ) as a separate equation using at the current iteration step 
a value for the conductive heat flux at the surface calculated from 
the quantities of the last iteration step. The resulting value of 
is prescribed as Dirichlet boundary condition, Eq. ( |33] |. for Eq. 
(l23T l. The temperatures calculated this way at each iteration 
step converge to the solution of Eq. (l23T l subject to Neumann 
boundary condition Eq. ( |32] |. This method worked without prob- 
lems. 

The time steps Af are chosen such that the relative change 
of the variables over At is about 3%. This is sufficiently small 
that a further reduction of the stepwidth does not significantly 
improve the accuracy of the numerical solution; a reduction of 
the admitted stepwidth by a factor of two results in our case in 
relative changes of the numerical values of the variables by a 
few times of 10"^. If the number of iteration steps becomes too 
big (e.g. > 20) with this choice, the stepwidth At is reduced by 
a factor of two until the number of iteration steps does no more 
exceed the limit. Since we use an implicit solution method, there 
is no limitation for the stepwidth from stability requirements. 

The initial model for the start-up of the solution method as- 
sumes a fixed temperature (= at initial time) within the body. 
An appropriate set of masses AM, is chosen that results (i) in a 
sufficiently fine grid at the surface to resolve the rapid tempera- 
ture variations at the surface and that (ii) is sufficiently fine for 
allowing to calculate the derivative dT/dMr at the centre with 
sufficient accuracy. For the initial model the porosities (p, and 
radii r,- for the set of mass-shells are calculated from hydrostatic 
equilibrium and the equation of state for cold pressing, as de- 
scribed in Sect. 14.21 

If a fixed temperature is to be prescribed at the outer 
boundary, this is technically achieved within the frame of our 
solution algorithm by letting Fgxt = crT^ in Eq. (l32l l. 

The solution method also allows to consider growing bodies 
by increasing the mass of the outermost shell according to some 
prescribed growth-law and splitting this shell into two shells at 
each instant where its mass exceeds some threshold value. This 
option in our code is not used, however, in the model calculations 
presented in this paper 

7.2. Some sample calculations 
7.2.1 . The model of Miyamoto et al. 

As a first test we calculate with the code a model using the same 
model parameters as in Miyamoto et aL (1981). The basic pa- 
rameters of the model are gi ven in Table [5]in t he col umn marked 
with MFT81. The model of iMivamoto et all (Il98lh is one for a 
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Fig. 8. Evolution of radial distribution of temperature for model 
PL (see Table |5] for its definition). Note the radial shrinking of 
the body by compaction of the initially porous material at about 
0.6 decay timescales of ^^Al. 



completely homogeneous body and does not consider the effects 
of porosity and the possibility of sintering. The data assumed for 
K and q correspond to average properties of L chondrites that, 
in fact, have low but non-vanishing porosities, scattering around 
about 4> = 10% (Yomogida & Matsui 1983). The true bulk den- 
sity and heat conductivity of completely consolidated chondritic 
material is higher (pb u ik = 3 .6 g cm"-', K - 4.27 W kg"' K"', see 
lYom ogida & Matsuil (Il983h . their Table 5). The model is calcu- 
lated by choosing as initial value (p - Q which guaranties that 
during the course of the calculation the porosity remains zero. 
Heating is only by decay of ^^Al. The result for the tempera- 
ture evolution in the centre of the body and a number of selected 
radii is show n in Fig. [7^. Th i s is al most identical with the result 
obtained by Mivamoto et al.l (Il98lh from the analytic solution of 
the heat conduction equation, i.e., our code reproduces this exact 
analytic result. 

7.2.2. Model of a porous body 

In Fig. [TJ? we show the results for the temperature evolution of 
a body having the same size and using a similar set of parame- 
ters, but now considering that the heat conductivity of the porous 
material, Eq. ( l30l l, is different from the value of heat conductiv- 
ity used by Mivamoto et al.l (1198 1'). and that the material from 
which the body forms is initially porous and consolidates by sin- 
tering. The parameters of the model are given in Table |5] in the 
column marked with PLO. 

It is assumed that the porosity of the surface layers at low 
pressures is 0srf - 0.6, co rresponding to the degree of com- 
paction found in lWeidling e t al. (2009) for powder material that 
was subject to numerous impacts. This is what one expects for 
the early formation time of asteroids where they grow by re- 
peated slow impacts of much smaller bodies. In deeper layers 
of the body where pressures are high due to self-gravity the ma- 
terial is compressed by isostatic pressing to higher densities up 
to a limiting value of <p ^ 0.4, see Sect. 14.21 The corresponding 
initial distribution of porosities in the interior is calculated as 
described in Sect 14. 21 For typical results see Fig.|2] This kind of 
compaction is a purely mechanical effect due to mutual rolling 
and gliding of the powder particles driven by an applied pres- 
sure which requires no elevated temperatures and acts therefore 
already in cold bodies (cold pressing). Starting from this ini- 
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Fig. 9. Evolution of radial distribution of filling factor D (relative 
density) for model PL (see Table |5] for its definition). Shown is 
the very initial phase of the evolution where the initially porous 
material is compacted by sintering. The resulting shrinking of 
the planetesimal size occurs at about 0.6 decay timescales of 
26 Al. 
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Fig. 10. Temperature evolution of test models for a porous bod- 
ies with modified values of parameters. Left panel with bound- 
ary temperature T\, = 150 K, right panel with T\, = 250 K. Upper 
panel with heat conduction coefficient K\, - 4W m"' K"' , lower 
panel with /Tb = 2 W m"' K"' . 



tial configuration the evolution of the model was calculated. The 
porosity dependence of the heat conduction is taken into account 
by means of approximation Eq. ( l30l l. The surface temperature 
was taken to be constant over time and equal to = 150 K. 

Figure [TJi shows that the peak temperature achieved in the 
centre of the body is l ower in this model than in the model of 
iMivamoto et al.l (|1981|) . This results from the higher heat con- 
ductivity in our model after complete sintering (Kb = 4 vs. 
K\, - 1). Contrary to this, the temperature in the outer lay- 
ers of the mode l increases more rapidly than in the model of 
IMivamoto et al.l (Il981i) because the high porosity outer layer 
acts as an insulating shield that prevents an efficient loss of heat 
to outer space. At a temperature of about 700-750 K (depending 
on pressure) sintering by hot pressing becomes active and the 
porosity rapidly approaches ^ = at higher temperatures. The 
temperature structure then becomes nearly isothermal in the inte- 
rior of the body and the temperature drops to the outer boundary 
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value within a layer of only a few km thickness. This is shown 
in detail in Fig. [8] 

Figure|9]shows the evolution of the porosity during the earli- 
est evolutionary phase. An outer shielding powder layer persists 
during the whole evolution of the body because cooling of the 
outer layers prevents the material to warm-up to the threshold 
temperature at about 750 K for compaction by hot pressing at 
low pressures. This behaviour i s com pletely analogous to what is 
found in Yomogida & Matsu|j 19841) and can be compared wit h 
results by .Sahiipal et al., (i2007h and iGupta & Sahiipall (l2010h . 
They considered gradual sintering in the temperature range of 
670-700 K by a smooth interpolation recipe, reducing the poros- 
ity from 55% to 0% by compaction of the body, and took into 
account respective changes in thermal diffusivity. 

Because the porosity approaches zero everywhere in the 
body where the temperature exceeds the threshold for sintering 
by hot pressing, the radius of the body shrinks significantly, typ- 
ically by 20% of its initial value. This occurs after about 0.6 
decay timescales of the main heat source, ^^Al, in this model, cf. 
Figs.[8]and|9] The size of the model, marked as PLO in Table|5] 
of 85 km corresponds to the final radius that the body would have 
after compaction to zero porosity (the initial radius before sin- 
tering is ^ 105 km). The final radius of the body almost exactly 
equals the hypothetical final radius at zero porosity, since the 
powder layer remaining at the surface is rather thin. Also the 
temperatures shown in Fig. [TJ? for a number of depth points be- 
low the surface correspond to that Lagrangean Mr-coordinate, 
which after compression to zero porosity would have the given 
value of the radius coordinate. Before the body shrinks these 
points have somewhat bigger depths below the surface. 

In Fig. [T]: we show the temperature evolution in a model 
with the same set of parameters as the previous model (see 
Table |5] model PL), that considers heating by decay of ^"Fe 
and long-lived radioactive nuclei as additional heat sources, us- 
ing an ^"Fe/''''Fe ratio at the upper limit of observed values (see 
Sect. 15.3b . The peak central temperature is about 30% higher 
than without this heat source. 



7.2.3. Mass-shells and time steps 

Since in the present models the mass of the body is constant, a 
fixed grid of mass-shells is used in the calculations. This grid 
was constructed as follows: Starting from an outer layer with 
some (small) thickness, the thickness of the layers from the out- 
side to the interior increases by a constant factor from one shell 
to next such that for some given number K of mass-shells the 
radius of the innermost sphere is 100 m; this fixes the width of 
the outermost layer The models of this paper are calculated with 
K - 300, in which case the outermost layer has a thickness of 
« 3 m. This choice of grid guarantees a sufficient high resolution 
of the thin outer region of rapid drop of temperature. An increase 
of K to 600 does not result in significant changes of the model 
characteristics; the relative change of central temperature, e.g., 
is * 10-4. 

In the model calculations used for fitting models to empiri- 
cal cooling histories of meteorites described in Sect. |8] the total 
number K of shells is increased io K - 1200. This is necessary 
if one requires that even in the region of steepest temperature 
decrease toward the surface the relative changes of temperature 
at some fixed mass-coordinate are at most of the order of a few 
times 10^4 if the number of grid points is doubled. 

The timesteps Af are chosen according to the strategy de- 
scribed in Sect. 17.11 The timesteps during the initial heat-up 
phase were typically a few thousand years. Once sintering com- 



mences, the step size decreases to about 100 years and varies 
around this value until complete sintering of the body (except 
for the outermost layers). Then Af increases continuously and 
during the final phase is limited to a maximum value of 1 Ma in 
order to obtain not too widely spaced steps for plotting purposes. 
The total number of timesteps required for a complete model run 
for an evolution period of 1 00 Ma is between 3 500 and 4 000. 



7.2.4. Dependence on parameters 

In order to get an impression how the model depends on some 
important parameters, we show in Fig. [TO] results for the same 
kind of model as model PL, but calculated with two different 
values of the surface temperature Ts and bulk heat conductiv- 
ity coefficient K^. The models in the left panel are calculated 
with a fixed surface temperature of 150 K, the models in the right 
panel with a fixed surface temperature of 250 K. These two val- 
ues encompass the possible values of the surface temperature for 
bodies in the region of the asteroid belt for both cases, if either 
the body is embedded in the accretion disc (cf. Fig. |2]i or the 
accretion disc has gone and the body is under the influence of 
the radiation of the young sun. There are only small differences 
between the run of the temperature at different depths, i.e., the 
temperature evolution below the immediate surface layer does 
not critically depend on the surface temperature, at least not for 
bodies with radii of the order of 100 km which are our main 
topic. Therefore we do not attempt in our further calculations 
to calculate Ts as precise as possible from Eq. (l32l i and simply 
assume a reasonable but fixed value for all times. 

The upper panel in Fig. [10] is calculated with a value for 
the heat conduction coefficient of K\j - 4Wm-' K"', the lower 
panel with - 2 Wm ' K The first value corresponds to 
what has been found as the average value for H and L chon- 
drites if measured values are extrapolated to zero porosity, see 
Sect. 15.41 As one can see from Fig.|6]there is significant scatter 
in the heat conduction coefficient (of presently unknown origin) 
and it is not clear whether the investigated sample of H and L 
chondrites are representative for the whole material of the parent 
body of the H or L chondrites. The value of - 4W m ' K"' 
corresponds to typical valu es of K for pure silicate minerals (cf. 
lYomogida & Matsuill98 3') and therefore probably represents the 
upper limit for the possible values of K^. Lower values, there- 
fore, may also be of interest for real planetesimals. Figure [TO] 
shows that the results of the model calculation depend signifi- 
cantly on the value of K^. Because it is presently not possible 
to determine from first principles for chondritic material, we 
consider K\j in our later model calculations to be a free param- 
eter (but, of course, restricted to the range of values found for 
chondrites). 



7.3. Maximum central temperatures 

The maximum central temperature that is reached during the 
course of the evolution of a planetesimal is an indicator of what 
kind of changes the material may undergo. If the central temper- 
atures exceed t he solidus t emperature of chondritic material of 
about 1 400 K jAgee et al.| [T995) partial melting occurs and the 
body starts to differentiate. If the temperature stays below the 
threshold temperature of about 700 K (at a; 0. 1 bar) for sinter- 
ing, the whole body retains its porous structure. The maximum 
central temperature Tc,max depends mainly 
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«form IMa] (form [Ma] tform [Ma] 

Fig. 11. Variation of maximum temperature in the centre of a planetesimal with radius R and instant of formation fform- Upper panel: 
Completely consolidated bodies with porosity = 0. Lower panel: Models of initially porous bodies with porosity (p^^ - 0.5 which 
consolidate in the interior during the course of their thermal evolution. The lines correspond to the indicated maximum central 
temperature. At temperatures above ^ 1 400 K partial melting of the mineral mixture is expected. Temperatures exceeding 1 500,K 
are left out for this reason. Shown are models for three different initial ''"Fe/^^Fe abundance ratios: Left panel: Without *"Fe. Middle 
panel: The optimum fit value of 4.1 x 10"^ for the H chondrite parent body, see Sect. [8] Right panel: Highest value of 1.6 x 10"* 
given in Table |5] 



1. on the radius, R, of the body, because this determines how 
efficiently heat can be removed from the central region by 
heat conduction, and 

2. on the formation time, fform, because this determines how 
much of the initial inventory of radioactive material already 
decayed before the body grew to its final size. 

Figure [TT] shows the dependence of max on R and fform for 
models of initially porous bodies and, for comparison, of bodies 
with pore-free material. Obviously the thermal evolution history 
of initially porous bodies is very different from history of equal 
sized compact bodies. Models are shown for three different as- 
sumptions on the abundance of ^''Fe as additional heat source 
besides ^*'A1. 

For porous bodies smaller than k 5 km radius the initial 
porosity is very high because they are even not compacted by 
cold pressing (cf. Fig.|2]). Because of their low initial heat con- 
ductivity even rather small bodies {R > 0.5 km) heat up at least 
to threshold temperature for sintering and become compacted in 
their interior, because the strongly insulating powder layer on 
their surface prevents their cooling. Completely compact bod- 
ies would reach central temperatures higher than 700 K only for 
radii exceeding a: 5 km because of much more efficient heat con- 
duction. 



For initially porous bodies bigger than » 50 km radius al- 
ready the initial porosity is low throughout almost all of the 
body because such bodies are already strongly compacted by 
cold pressing (cf. Fig.|2| and the remaining porosity rapidly dis- 
appears by sintering. Their thermal evolution is essentially iden- 
tical to that of completely compact bodies, except, of course, 
in the layers close to the surface that retain part of their initial 
porosity. 

Porous bodies with radii between ^ 5 km and » 20 km are al- 
ready significantly compacted in their central part by cold press- 
ing (cf. Fig.|2]l but have initially an extended low-porosity outer 
mantle. Porous bodies with radii between =s 20 km and ^ 50 km 
also are already compacted throughout the body by cold press- 
ing (cf. Fig.|2]), except for the outer » 10% of their radius. They 
show the most complex dependence of Tc,max on R and formation 
time. 

Temperatures above rc,maxl500K are not considered be- 
cause we consider in this paper parent bodies of undifferentiated 
meteorites. A t a temperature o f T^max <: 1400K the silicates 
start melting dAgee et aLlll995h and differentiation is unavoid- 
able. 
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Fig. 12. Optimum-fit model for the cooling 
history of the parent body of H chondrites. 
Abscissa is time elapsed since formation CAIs. 
Evolution of temperature at a number of depths 
below the surface is shown. The upper con- 
tours of shaded areas correspond (from bottom 
to top) to depth's of 0.32 km, 2.3 km, 7.8 km, 
1 1 km, and the uppermost contour to the centre. 
The rectangular boxes and circles correspond 
to the empirical data for H6 and H5 chondrites, 
respectively, given in Table|6] Crosses are error 
bars. The dashed lines correspond to the tem- 
perature evolution at depth's of 8.9 km (lower 
line) and 36 km (upper line) below the surface. 
They represent our best fit to the empirical data. 



Table 6. Key data for cooling history of selected H chondrites 



Quantity 


Kernouve Richardton 




type 


H6 H5 






Hf-W' (metal-silicate) 




radiometric age 


4559±1 4563±1 


Ma 


temperature 


1150±75 1100±75 


K 




U-Pb-Pb' (phosphates) 




radiometric age 


4522+1 455 1±1 


Ma 


temperature 


720±50 720±50 


K 




Ar-Ar^-' (feldspar) 




radiometric age 


4499±6 4525±11 


Ma 


temperature 


550±20 550±20 


K 




Pu-fission tracks (pyroxene-merrilite) 


calculated age'' 


4499±6 4525±11 


Ma 


cooling rate 


2.6±0.5 2.9±0.5 


K/Ma 


calculated age'* 


4438±10 4469±14 


Ma 


time interval'-^ 


61±8 56±9 


Ma 




metallographic 




cooling rate'-' 


10 20 


K/Ma 



Notes: 

' Time -interval for Pu-fission track cooling rate from 550-390 K, metallographic cooling 
rate 800-600 K 

^ Recalculated for miscalibration of K decay constant (explanation see text and 

|Trieloffetalj200l 

^ Data f rom Hf-W: lKleine eTaP Oool U-Pb-Pb: JGopel et aU1994l) . Ar-Ar: iTrieloff et al] 

I2003D , metallographic cooling rates:|Ta vlor et alj|l987|) 
^ Calculated age at 390 K from time interval between Pu-fission track (merrilite, 390 K) and 

Pu-fission tracks at 550 K (pyroxene) compared to Ar-Ar feldspar age at 550 K 



8. Thermal history of the H chondrite parent body 

8.1. Empirical data for cooling history 

Most meteorites reach the Earth as cm- or m-sized rocks, as they 
are the result of repeated impact fragmentation of the initially 
much larger parent bodies. As such events can be highly ener- 
getic, they change both chemistry, structure, and also disturb the 
isotopic memory (i.e., the age information). Hence, information 



Table 7. Properties of optimum fit model 



Quantity 




value 


unit 


Radius 


R 


100 


km 


formation time 


^form 


2.3 


Ma 


heat conductivity 




4.0 


Wm' K' 


surface temperature 




300 


K 


60pg|56pg aijuudance ratio 




4.1 X 10-'' 




initial porosity 




50% 




Chondrite type 


depth 


range 






from 


to 




powder layer 


0.00 


0.29 


km 


H3 


0.29 


2.3 


km 


H4 


2.3 


6.7 


km 


H5 


6.7 


10.8 


km 


H6 


10.8 


100 


km 


burial depth 


Kernouve 




36 


km 


Richardton 




8.3 


km 



of cooling histories extending back to the origin of the solar sys- 
tem must be restricted to meteorites that 

1 . show extraordinarily low mineralogical and structural char- 
acteristics of shock metamorphism induced by asteroid col- 
lisions, 

2. were dated with high precision, and 

3. were dated simultaneously by a set of various high and low 
temperature chronometers tracing the cooling history from 
high temperatures (> 1 000 K, e.g., Hf-W) down to interme- 
diate (e.g. U-Pb-Pb or K-Ar ) and very low temperatures (if 
possible < 400 K, e.g. ^'^"^Pu fission tracks). 

Such high quality data on un-shocked chondrites are quite lim- 
ited, in spite of the high number of meteorites available in ter- 
restrial collections. While in the case of L cho ndrites, a major 
impact 470 Ma ago dKorochantseva et al]|2007) seems to have 
deeply erased the primordial low temperature cooling history, H 
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chondrites seem more promising. For a number of seven - notice- 
ably un-shocked - H chondrites, complete high precision Hf-W, 
U-Pb-Pb, Ar-Ar and ^"^'^Pu fission track data along with metallo- 
graphic cooling rate data exist. Table |6] shows such data for the 
H6 chondrite Kernouve and the H5 chondrite Richardton, which 
we can use for a preliminary sample calculation. 

8.2. Fit of selected H chondrites 

A "fit" to the data in Table |6] is shown in Fig. [12] The chrono- 
logical data for these two meteorites best fit cooling curves 
in an asteroid at 36 and 8.9 km depth. The properties of the 
parent body in Table |7] were obtained according to the fol- 
lowing fit-procedure: The initial abundance of ^^Al is roughly 
constrained by the ^^A1/-^A1 ratio of ordinary chondrite chon- 
drules, which is an upper limit shortly before parent body for- 
mation. Chondrule measurements indicate an ^^Al abundance 
corresponding roughly to 2-3 Ma formation time after CAIs. 
Similarly, the ^"Fe abundance is constrained by primitive type 3 
ordinary chondrites (sulphides, see iTachibana & HussI |2003|) . 
Furthermore, abundances of ^^Al and ^''Fe (or, in other words, 
the formation time of the parent body) must be such that suffi- 
ciently high maximum metamorphic temperatures in the aster- 
oid are achieved to allow strong type 6 metamorphism, but also 
not too high to prevent partial melting, metal silicate separation 
and differentiation. Heat conductivity and radius mainly influ- 
ence the total duration of the cooling time of the parent body. 
We arbitrarily chose 100 km, although a smaller asteroid would 
also allow extended cooling as observed for Kernouve, which 
in this model needs only a burial depth of about 36 km. The 
boundaries separating H6 from H5, H4 and H3 material are rel- 
atively shallow in the model, due to the insulating porous thin 
outer layer 

8.3. Discussion 

The most prominent feature in our new model is the possibil- 
ity of relatively small parent planetesimals with significant heat 
retention. In the H chondrite parent body model, this shows up 
in relatively thin layers of less heated or metamorphosed ma- 
terial. Moreover, the relatively fast cooling required to achieve 
temperatures below Hf-W and U-Pb-Pb closure for the H5 chon- 
drite Richardton (within 3 and 13 Ma, respectively) sets an upper 
limit to the contribution of decay heat of ^"Fe (roughly about 20- 
30%). For example, if ^°Fe contributed all decay heat of ordinary 
chondrite parent body metamorphism, suflicient fast cooling of 
H5 metamorphosed material would not have been possible, as 
the half -live of ^"Fe is 2.6 Ma (a new revised value) versus 
0.72 Ma for ^^Al, implying significantly longer heat production 
and implicit slower cooling in the first few Ma. This result is in 
line with the initial ^ °Fe concentration found in primitive type 3 
ordinary chondrites (ITachibana & Husi l2003l). lower than initia l 
values previously obtained for CAIs (Birck & Lugmaiij[l988h . 
and supports the view that ^"Fe was likely not distributed homo- 
geneously in the early solar system. A more detailed H chondrite 
parent body modeling will be presented elsewhere. 



9. Concluding remarks 

We constructed in this paper a model for the thermal evolu- 
tion of the parent body of H chondrites. The model consid- 
ers compaction by cold pressing and sintering by 'hot press- 
ing'. The heat conductivity of the porous material was de- 



termined by combining new data obtained by iKrause et al.l 
(1201 lahf or high porosity material with data for porous chon- 
drites (jYomogida & Matsui 1983). A model was fitted to data 
on the cooling history for two H chondrites, Kernouve (H6) and 
Richardton (H5), for which particular accurate data are avail- 
able. It was shown that it is possible to find a consistent fit for the 
parent body size and formation time that reproduces with good 
accuracy the empirically determined cooling history of both, H5 
and H6 chondrites. 

For obtaining our consistent fit, it was necessary to include 
(besides radius of the body and formation time) also the abun- 
dance of ''"Fe into the fitting procedure. A value of ^''Fe/''^Fe was 
determined, that is within the range of values reported in the lit- 
erature for different meteorites. No other arbitrary fit parameters 
are required; all other properties of the model are fixed by the 
physics of the problem. 

The new model predicts rather shallow outer layers for petro- 
logic types 3 to 5 because of the strong blanketing effect of an 
outer powder layer, that escapes sintering. In so far the model 
deviates considerable from previous models that are based on 
the analytic model of (Mivamoto et al. 1981). Other properties 
of the model are similar to older models; in particular radius and 
formation time are not substantially different. 

The present model, though relaxing some earlier shortcom- 
ings, still has a number of shortcomings. The most stringent is 
the instantaneous formation hypothesis, that needs to be relaxed 
because the formation time of the body by run-away growth (of 
the order of 10^ a) is shorter than the decay time of ^^Al, the 
main heating source, but probably not short enough for being 
completely negligible. The second severe shortcoming is that 
heat conductivity of porous media can not yet be treated from 
first principles on. This is not possible with presently available 
computing resources, but this may change in the near future. 
Other shortcomings are a rather simplistic treatment of sinter- 
ing and of the outer boundary temperature. We modelled in this 
paper sintering with the sam e kind of theoretical description as 
lYomogida & Matsuil (1 19841) in order that our results can he com- 
pared with that model. This treatment of sintering should be re- 
placed by more elaborated modern theories of hot pressing. 
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